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ABSTRACT 

We attempt to measure the density of dark matter in the two Dwarf Spheroidal Galax- 
ies Fornax and Sculptor using a new method which employs Jeans equations based 
on both the second and fourth moment of the Collisionless Boltzmann Equation (i.e. 
variance and kurtosis of line of sight stellar velocities). Unlike previous related efforts, 
we allow the anisotropy of the radial and tangential second and fourth order moments 
to vary independently of each other. We apply the method to simulated data and es- 
tablish that to some degree it appears to be able to break the degeneracies associated 
with second order only Jeans analyses. When we apply the technique to real data, we 
see no huge improvement in our understanding of the inner density of Fornax, which 
can still be fit by either a quite cuspy or cored density profile. For Sculptor however 
we find that the technique suggests that the data is incompatible with a steep profile 
and a cored profile seems more consistent. As well as presenting these new results and 
comparing them to other estimates in the literature, we try to understand why the 
technique is more effective in one system than the other. 
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1 INTRODUCTION 

It is well documented that the ACDM model of cosmol- 
ogy is remarkably successful at reproducing the observed 
large scale structure of the universe (Tcgmark ct al. 2004;) 
in dissi pationless N-body si mulations of cold dark mat- 
ter (e.g ISpringel et all [200i). At galactic scales however 
there are a num ber of discrepancies between the simula- 
tion pred ictions (^ Navarro et al.l Il996l : ISpringel et al.l l2008l : 
iKlvpin et al. 2011i) and ob servations inc l uding (but not lim- 
ited to) the abundance l|Klvpin et al.l [l999l: iMoore et al.l 
19991) and phas e space correlation (|Kroupa et al.l l2005l : 



Conn et al.ll20ld ) of satelli te galaxies, the absence of h ighly 



luminous satellite galaxies (|Bovla n-Kolchin et al."201ll) and 
the observation of seemingly cored de nsity profiles at the 



centre of low surface brightnes s spiral llGentile et al 



and dwarf spheroidal gala xies l|Walker fc Penarrubia 



2004) 



2011 



lAmorisco fc Evang l2012al ) that is at odds with the more 
cusped density profiles ubiquitous in the halos of cold 
dark matter simulations. The interaction of dark matter 
with baryons, absent in most simulations, could ameliorate 
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llTonini et al.ll2006l : iGovernato et al.ll2012l : iRead fc Gilmord 
I2OO5I ) any of these differences to a greater or lesser ex- 
tent and is an area of active study. Work has also been 
undertaken to investigate the possibility that the Milky 
Way is a statistical out lier l|Strigari fc Wechslerl I2OI2I : 
iBovlan-Kolchin et al.ll2010l '). 

The interesting anomalies have also given rise to a 
number of alternative cosmological models. Warm dark 
matter is one such contender with a free-streaming 
length that pr eserves the successful large-scale formation 
llBovarskv et al..i,2009.') wh ilst potenti ally reducing the ab un- 
dance l|Colm et al.l I2OO0I ) and size l|Lovell et al'l |2012| ) of 
satellites in line with observation. It has also been suggested 
(Zavala et al. 2013) th at relaxing the assumpti on of colli- 
sionless dark matter (jSpergel fc Steinhardt|[200oh could gen- 
erate extended density cores. 

To minimise baryonic effects, the high mass- luminosity 
ratio observed in dwarf spheroidal galaxies makes them a rel- 
atively 'clean' target to test predictions from ACDM simula- 
tions and herein we use satellites of the Milky Way to try and 
tackle the cusp vs core issue. This property also makes them 
good candidates for indirect detection sea rches of dark mat- 
ter su ch as that with the Fermi telescope l|Ackermann et abl 
l201lh and a precise measurement of the density profile is 



© 0000 RAS 



2 T. Richardson & M. Fairhairn 



key to reducing the large uncertainties in the the expected 
flux of incident annihilation products such as gamma rays. 

Several methods have been proposed to infer the density 
of dwarf spheroidal galaxies from a limited sample of pro- 
jected radii and velocity measurements. A key difficulty that 
the majority of these methods face is that to infer the gravi- 
tational potential one must contend with the six dimensional 
phase space distribution function /(r, v) as an unknown nui- 
sance parameter which is further obscured from the observed 
velocity data by projection along the line of sight. 

P erhaps the simplest metho d is to use the Jeans equa- 
tions l|Binnev fc Tremaindliooj ) that relate the underlying 
gravity potential to the more manageable moments of the 
phase space distribution function. As partial solutions to 
the collisionless Boltzmann equation however then with a 
truncated series one cannot guarantee that all solutions cor- 
respond to physica l distribution functions and must contend 
with degeneracies (iMerrittlligSTi ) upon projection along the 
line of sight. The degeneracies associated with the tradi- 
tional second order Jeans analysis are more pronounced for 
a flat stellar velocity dispersion profile. Being an analytic 
method it is also necessary to parametrise the anisotropy 
and density that for mathematical simplicity alone have led 
to constraining and often unphysical additional assumptions 
such as Gaussianity and constant anisotropy in the litera- 
ture. 

To alleviate issues with the Jeans equations raised 
above, n umerical procedures such as the orbit-superposition 
method (|Schwarzschildlll979l ) utilise the Jeans theorem to 
ensure physical solutions to the collisionless Boltzmann 
equation and without any assumed form for the anisotropy. 
Such methods have not yet however been able to dis- 
tinguish cusps from cores in classical dwarf spheroidals 
l|Breddels et al.ll2012l : Ijardel fc Gebhardtll2012l). By c onsid - 



ering the simpler analytic Jeans analysis, Wolf et al] (|2010l ) 
discovered that the mass at the half-light radius is invari- 
ant under a general choice of anisotropy. With the addi- 
tional discovery of multiple stellar populations in Fornax 
and Scul ptor dwarf spheroidals this was used to devise a 
method (jWalker fc Penarrubial I2OIII ) that places perhaps 
the strongest constraints in favour of cored density profiles 
in dwarf spheroidals. 

To test the assumptions mentioned in the paragraph 
above that underpin this result and to alleviate the degen- 
eracy that blights a dispersion only method we were moti- 
vated to extend the truncated Jeans analysis to fourth or- 
der. In the context o f dwarf spheroidals this w as first pro- 
posed in iLokad (I2OO2I ) and was applied to Draco iLokas et al] 
l|2005l ) but was limited to a constant anisotropy parame- 
ter and phase space distribution functions with correlated 
second and fourth moments where tangentially biased vari- 
ances necessa rily implied tangentially biase d kurtosis. In 
previous work I Richardson fc FairbairnI l|2012t) we provided a 
framework that, by introducing anisotropy parameters anal- 
ogous to Binney's at second order, imposed no assumptions 
on the form or correlations of anisotropy parameters. This 
extension facilitates the most generalised representation of 
anisotropy yet used in an analytic Jeans analysis that com- 
plements numeric approaches. 

A brief summary of the fourth order Jeans analysis and 
the statistical likelihood analysis used to fit this to the ve- 
locity data is provided in Sections 2 and 3 and we refer 



the reader to iRichardson fc" FairbairnI l|2012l ) for a detailed 
account. In Section 4 we display the results of the analy- 
sis on data sets from Fornax and Sculptor dwarf spheroidal 
galaxies. Final we discuss the implications of the results in 
comparison to other works and make concluding remarks. 



2 SOLUTIONS TO THE JEANS EQUATIONS 
AT FOURTH ORDER 

One particular means to infer the DM dominated density 
of a dwarf spheroidal galaxy from a set of stellar posi- 
tions and velocities is th e collsionless Boltzmann equation 
(jMerrifield fc Ke"n3ll990l ) which for equilibrated spherically 
symmetric systems is 



dt dr \ r dr j dvr 



df 



r 
1 
r 
0. 



cote — VrVe) 
[v^Vr + V^Vg cot ( 



a/ 

dve 



(1) 



df 



where f{r,w) is the four dimensional distribution function 
of stellar radii and velocity components in spherical coordi- 
nates Vr,ve,v^ and ^{r) is the spherically symmetric grav- 
itational potential that we assume to be dominated by the 
dark mass component. Equation ^ alone is not easy to 
work with but by multiplying through with a suitable factor 
and integrating over all velocities (^Merrifield fc Kent .l990h 
one obtains the Jeans equations that stipulate conditions on 
the moments of the distribution for equilibrated systems. 



Oi 2i 2k / 2i 2j 2k ri \ 

uv^'-Vg-'v^'^ = / Vr Vg-'v^ f(r,v)d v. 
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of which v{r), the zeroth moment that effectively normalises 
the distribution function is the observed lo cal stellar density 
and the others may be written l|Deionghe|[l986.1 as 



"r 



2(j + k) 



(3) 



where B{x, y) is the Beta function and vt = y + «^ is the 
magnitude of the two dimensional tangential velocity. Due 
to the assumption of spherical symmetry net rotations and 
all other odd moments vanish. The two moments at second 
order, i.e the variances = «^,o"t = axe related to the 
underlying gravitation al potential by the w ell-known second 
order Jeans equation (|Binnev fc Tremaind [2008 1. 

d{vaf.) 



dr 



(4) 



At fourth order t here are two unique Jeans equations 
((Merrifield fc Kentlflgga ) that relate the three fourth order 
moments, 



d{vv^ 
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With n equations and n + 1 moments at the 2n"' order then 
even if one has knowledge of the stellar and DM densities 
there still persists a degeneracy of solutions to the Jeans 
equations at each order. Upon projection along the line of 
sight this unknowable choice gives rise to the degeneracy 
problem. To represent the full set of degenerate solutions 
the parameter, 



/?(r) = 1 - 



(7) 



2a2(r)' 

which measures the second orde r deviation of the system 
from isotropyQ, was introduced (|Binnev fc Tremainell200^ ) 
to fully constrain Q and its specification defines the solu- 
tion. We replicate this idea at fourth order by introducing 
an analog to Binney's anisotropy parameter 



(8) 



that represents the deviation of the fourth order moments 
from the isotropic system v$ = f^'fr ? and therefore naively 
(see iRichardson fc FairbairnI |2012| . for details) determines 
the difference between the kurtosis of the radial and angular 
velocity distributions. All solutions to the Jeans equations 
at fourth order may then be represented by the ordinary 
differential equations, 



d{. 



2/3 2^ 
ar r ar 



0. 



dr 



H uv^ + iua^— — 0, 



dr 



(9) 



(10) 



from which one can deduce the remaining moments with 



2(1 - /3)cr?, v^v^ = 1(1 - ^')«r and 



(15) 

where R is the observed, projected radius perpendicular to 
the line of sight. The function g{j3',r,R) common to the 
separable augmented density system is 



g{P',r,R) = l 



- /3') i?^ 



-4^3^ (16) 



and w e note the general Jeans solution extends the work of 
iLokad (|2002l ) for which /? = /3 = constant and its gener- 
alisation to arbitrary ani sotropy, t he separable augmented 
density model (see for e.g. lAn|[201lf ) defined by /3'(r) = /3(r), 
with an additional term. These projected moments provide 
a means of comparison with observable data alongside the 
surface density 



E(i?) = 2 



R^ 



I dr. 



(17) 



Now we are in possession with expressions for the second 
and fourth moments of the line of sight velocity dispersion 
for a given underlying dark matter density, we need to see 
how well we can constrain the unknown dark matter density 
using data. 



3 LIKELIHOOD ANALYSIS 

In this section we briefly summarise the Jeans/Monte- 
Carlo Markov chain (MCMC ) methodology presented in 
IRichardson fc FairbairnI (|2012l ') by which the posterior dis- 
tributions for a set of anisotropy and density parameters p 
are generated from a data set d of stellar line-of-sight radii 
and velocities via likelihood function £{d\p). 



U,^ + 2(/3-/3)ra.^-. 



(11) 

If one assumes the anisotro py to be co nstant with radius 
then an integral form exists ()Lokasll2002l ) that simplifies the 
calculation, 

iya^.{l3 = const) = r"^'' H i^^u^dr (12) 



v^[j3' — const) = 3r' 



-2/3' 



' 2/3' 2d$ , 

r va^ — dr 
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although in this work we will be considering situations where 
both anisotropy parameters vary with radius. Once we are 
in possession of all the necessary moments it is possible to 
project along the line of sight and calculate the moments of 
the line-of-sight velocity distribution, 



Eap(i?) = 2 



R^ 



zdr 



(14) 



/OO 
(3(/^')^+l^(/3'-/?)'^r^^ 



v{r)r 



i?2 



zdr 



^ i.e the difference in widths of the radial and angular velocity 
distributions 



3.1 Anisotropy and Density Parameterisations 

With the ultimate aim of quantifying the uncertainty in a 
measurement of a dwarf spheroidal's density profile a para- 
metric form for the anisotropy and density is required for an 
MCMC analysis. In particular we are motivated to consider 
realistic density profiles which can exhibit either extended 
cores or the cuspy centres predicted by dissipationless DM 
only simulations. To this end the Einasto parametrisation 
(Einasto & Haud 1989) of the density and its corresponding 
logarithmic density slope are 



Pdm(j") = p-2 exp { 



r 

r-2 



7(r) 



din p 
din r 



r 
r-2 



(18) 



(19) 



where r_2 is a scale radius that indicates where the loga- 
rithmic density slope y{r) = —2, p_2 is the scale density at 
this radius and a is a constant shape parameter that deter- 
mines the rate at which the density slope deviates from -2 
at the scale radius. Though all Einasto profiles are cored at 
the galactic centre 7(r — >■ 0) = 0, variation in the shape pa- 
rameter a can mimic cusped profiles at all but vanishingly 
small radii. If the shape parameter is small a < 1 then 7 
varies slowly from the scale radius r-2 to the centre of the 
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galaxy, extending the steeper gradient and wiping out the 
core. An estimate of the core size can be defined as, 

f 



lo: 



(20) 



which indicates the radius at which from zero the slope falls 
to 7 = —0.2 and we note that increasing a shifts the core 
radi us towards r-2. To mod el the anisotropy parameters we 
use (IBaes fc van Hesell2007l ). 



/?(r) = (/?oo - /?o) 



+ /3o 



(21) 



which generically describes anisotropy with a quadratic 
transition about rp that asymptotes at /3o for r = and /3oo 
for r — ^ 00. We will use the same functional form for the ra- 
dial dependence of /9'(r), with all parameters corresponding 
to /3' denoted with a prime. At second order the parameter 
set is thus p2 = {/3o, /3oo, J'/?, P-2, r_2, a} which is extended 
at fourth order by the additional fourth order anisotropy 
parameters, pi = {Pq, P^.rp, P'q, P'^,r'p, p-2,r-2,a}. 

The local density v{r) is modelled by a Plummer profile. 



u{r) oc 



Rl 



/2 



(22) 



where R1/2 is the Plummer radius which is equivalent to 
the projected half-light radius that encompasses half of the 
stars on the circular projected surface. This parameter is 
taken from observational data and we are motivated to use 
this simple approximation to reduce the already considerable 
number of free parameters in the analysis. 

We note that this is the most rigid part of the analy- 
sis and that imposing a central core for the stellar density 
can have a pronounced impact on the dispers ion at small 
radii t hat limits the range of physical solutions l| Evans et al.l 
Hooi)). Assuming that the dark matter and tracer den- 
sity have central slopes shallower than the isothermal cusp 



p(r) 



at the galactic centre and also that the 



anisotropy parameter asymptotes to a flat constant value 
/3o as modelled above, then the anisotropy term dominates 
the mass term such that. 



7* 



2A) 2 



(23) 



dr 

where 7* = —d\nv / d\nr is the logarithmic density slope 
of the tracer stars. With these assumptions then, that en- 
compass cored and NFW dark matter density profiles with 
logarithmic slopes of 7 = and 7=1 respectively the dis- 
persion is, 

lima.^ocr^*-^''". (24) 

r— fO 

The seemingly benign assumption of a cored Plummer pro- 
file and spherical symmetry therefore demands that the dis- 
persions of tangentially biased systems /So < vanish and 
whilst for radially biased systems /3o > they diverge at the 
centre. One must therefore take care when interpreting re- 
sults near the galactic centre and clearly the assumption of 
constant anisotropy at all radii coupled with a fixed tracer 
density is liable to these assumptions. Therefore though we 
adopt a Plummer profile with a fixed density slope the free- 
dom of /3o in the generalised parametrisation of jS effectively 
absorbs 7*, minimising the impact of the plummer profile 
assumption to some extent . 



3.2 Sample Moments and Likelihood Function 

With the Jeans equations defined in Section 2 it is possible 
to obtain the dispersion and kurtosis Kp of the line of sight 
velocity distribution corresponding to a set p of anisotropy 
and density parameters. We then fit these to the data set 
d of Af* stellar line of sight velocities {vi} and positions 
{Ri} with a suitable likelihood C{d\p) which denotes the 
probability that d is drawn from a phase space distribution 
function with moments <Jp{R) and Kp{R). Whilst one max- 
imises the information in the data set by defining the total 
likelihood as the product of individual tracer probabilities, 
£.{d\p) = ^{'^iWp{Ri)t f^p{Ri))j it is difficult to explic- 
itly encompass the highly significant statistical fluctuations 
from limited sampling into the probability distributions J-. 
As such we split the data into Nr radial bins of equal stellar 
content A'' = N^/N,. and use the sample moment estimators. 



1 2 



N 



1 

^= A^ 



("i - 

(52)2 



(25) 



(26) 



where fi is the mean velocity of the bin, to represent the ve- 
locity data. In this way, wherein a compromise is made be- 
tween statistical precision and spatial resolution, the proba- 
bility distribution Sx{x\ap{Rhin), «p(i?bin)) of the estimator 
X has a dependence on the sample size that can be calcu- 
lated numerically by simply taking many bootstrap sam- 
ples of size N from a suitable choice of parent distribution 
J-{v\ap{R]^in), K,p{Rbin)) and adding Gaussian noise to rep- 
resent the experimental errors. The distribution of kurtosis 
estimators was found to be both significantly biased and 
skewed when realistic samples of A < 400 velocity mea- 
surements were drawn from parent distributions for which 
the variance and kurtosis are known (we use both Pearson 
and Gaussian superposition as our non-Gaussian profiles). 
To better approxi mate Gaussianity the kurtosis estimator 
was transformed (|Lokas et al ] |2005l ) to K.' = (In k) ' but 
a significant skewness persisted for leptokurtic parent dis- 
tributions. For this reason a form for the likelihood is 
inappropriate and we instead took of order lO" bootstrap 
samples to generate a smootl|f|, fully numeric probability 
distribution by binning the sample kurtosis data. 

For simplicity the variance and kurtosis estimators were 
considered independent in our analysis and we found that 
the marginalised kurtosis distribution 5^/ is almost com- 
pletely independent of the input variance, so we assume it 
is entirely independent. With these assumptions we employ 
the following likelihood for the joint analysis of dispersion 



^ Maximum likelihood estimators were also considered but whilst 
they reduced bias and statistical errors the sampling distributions 
also displayed a significant skew. Maximising the 2D parameter 
space proved too computationally expensive to generate sufficient 
samples for a smooth numerical distribution. 

Though upon measurement the estimators are in fact weakly 
correlated with Pearson coefficient p 0.2 
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and kurtosis, 



£(d|p) = Y[S^{^^\aliR,),K;{R,)) X S.iKlK^iR,)). (27) 



and for comparison we use the dispersion-only likelihood, 



JZ^(d\p) = Y[Sa{^f\crl{Ri),Kp = 3). 



(28) 



for which the parent line-of-sight velocity distributions 
are assumed to be Gaussian. In practice we use the Pear- 
son family of distributions to generate the sampling dis- 
tributions >S for the joi n t ana lysis and refer the reader to 
iRichardson fc FairbairnI (j2012l ) for a visual and mathemati- 
cal guide to these distributions and a more detailed outline 
of the numerical procedure. 



3.3 Monte-Carlo Markov Chain Methodology 

An efficient means of deriving posterior distributions for a 
highly dimensional data set p is to employ a random walk 
Monte-Carlo Markov Chain (MCMC) that moves through 
the parameter space by proposing new positions p' in the 
chain according to the current location p and accepting 
these positions accordin g to the ratio of likeli hoods with 
the Metropolis-Hastings (|Metropolis et al]ll953l ') algorithm. 
An acceptance ratio of 0.2 to 0.4 strikes a balance between 
exploration and precision suitable for higher dimensional pa- 
rameter spaces. We adopt a multivariate Gaussian proposal 
density Q{p'\p) wherein the covariance matrix is updated pe- 
riodically from previous entries in the chain to optimise mix- 
ing of correlated parameters. A diagonal covariance matrix 
is used as an initial estimate where entries are determined 
by fixing all but one parameters and adjusting the effective 
univariate Gaussians width until an acceptance of « 50% 
is achieved. As the MCMC is maximally efficient when the 
proposal density matche s the posterior dist ributions the fol- 
lowing transformations jCharbonnier et al.. 201 1) and prior 
ranges are imposed on the anisotropy and density parame- 
ters, 



logio[l-^] 
log 10 B 

logio P-2 



log 



loe 



[-1,1] 
[1,4] 
[-8,5] 
[1,6] 

hl.3,1.3] 



(29) 



where A = {Po, Poo, P'o, P'^o} and B = {r,3,r^}. We run four 
MCMC chains for each analysis simultaneously from differ- 
ent starting positions and stop when the four MCMC chains 
satisfy the Gel mans- Rubins convergenc e criterion R < 1.03. 
We also check (|Charbonnier et al.ll201lh that the parameter 
variances across the chains are smaller than 1% of the mean 
of these variances. Inspecting the autocorrelation function 
of each chain we thin them out accordingly after removing 
the initial burn-in period. 

From the processed chains we are able to construct 
probability distributions and thus confidence intervals for 
a number of interesting physical quantities such as the 
enclosed mass and the logarithmic mass slope F = 
dlnM/dlnr. To do this a chain is created for each in a 



range of radii by evaluating these quantities from the ele- 
ments in the parameter chain. Assuming flat priors on the 
derived quantities the median values then correspond to best 
estimates and one may use fractiles to determine the cen- 
tral confldence intervals, i.e if 95% of F(200pc) values in the 
chain lie between F„ and Td then these become the 95% 
confldence intervals at this radius. 



3.4 Test Data 

Though the statistical method presented above has 
been tested ( {Rich ardson & Fairbairn 2012) for constant 
anisotropy, the addition of four parameters to the gener- 
alised case adds an additional strain to the MCMC analy- 
sis that warra nts further testing. Simu l ated d ata was gen- 
erated as per IRichardson fc FairbairnI (|201 2^ from a par- 
ent distributior0 with dispersion and kurtosis profiles cor- 
responding to parameter set pin which we aim to recover 
with the MCMC. As expected the additional parameters in- 
creased the number of iterations typically required to satisfy 
the convergence criterion and chains with between and 10^ 
and 10^ entries were thinned to one in every hundred for 
satisfactory autocorrelation values and independence. This 
could be improved in future analysis by adopting a more 
fiexible proposal density that can identify the different cor- 
relations between parameters in different parts of the pa- 
rameter space. 

Fig. [T] shows the moment estimators of the simulated 
data and the moment profiles corresponding to pin ~ {Po = 
0.25, /3oo = -0.2, = O.Skpc, = 0.0, pi^ = 0.2, rjj = 
0.2kpc,p_2 = O.O5M0pc"^,r_2 = Ikpc, a = 3.5} in the 
Jeans equations that is used to generate it. The error bars 
and central values show the confldence intervals and bias 
from the sampling distributions 5 (pin). Figs. [2] and [3] show 
that the MCMC was able to reproduce the input parameters 
Pin within 95% confidence intervals for a mock Fornax data 
set of two thousand stars, a plummer profile with a projected 
half-light radius of R1/2 ~ 575pc and normally distributed 
experimental errors of 22% relative to the intrinsic disper- 
sion equivalent to a few kms~^. The input parameters were 
chosen such that the density had an extended core and the 
dispersion measurements are approximately constant at all 
radii as shown in Fig. [1] with the curved features masked by 
the natural scatter. Figs. [2] and show that the dispersion- 
only fit is not able to distinguish between radial and tangen- 
tial anisotropy, nor is it able to tell the difference between 
NFW-like and cored halos. 

With the addition of the kurtosis measurement this de- 
generacy is broken with the dotted NFW profiles disfavoured 
to high significance and the cored nature of the simulated 
halo being detected. This shows that in some circumstances 
it is possible with generalised as well as constant anisotropy 
to break the degeneracy in the traditional analysis. This re- 
sult is of course aided by the chosen set of parameters with a 
kurtosis profile that has a distinct leptokurtic feature around 
300pc that disfavours platykurtic («:' < ln(3)^''^'' ^ 1.01) 



* Data was generated from th e Gaussian Superposition family of 
distributions in appendix B of iRichardson &: FairbairnI l l 20121) so 
as to again test the assumption of Pearson generated sampling 
distributions 
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Figure 1. Mock Fornax simulated data: Solid black curves indi- 
cate the the line of sight dispersion II14I I and kurtosis JTSj used 
to generate the mock data set pi^. Green triangles show the dis- 
persion l(25j and kurtosis of the mock data set. Red circles and 
error bars show the median and 95% confidence interval for the 
reconstruction of the inputc parameters following the analysis in 
the text. 

points in that region. With a perfect assessment of the tracer 
density u{r) and exact parametric forms for the anisotropy 
and density the simulation is always likely to be optimistic 
in its assessment of the confidence intervals but serves as a 
proof of concept and a crucial test of the MCMC's reliability. 



4 APPLICATION TO FORNAX AND 
SCULPTOR DWARF SPHEROIDALS 

4.1 Data 

Stellar position and velocity data for Fornax and Sculp - 
tor is taken from the Magellan survey (| Walker et al.l [20091 ) 
and after using the spectral metallicity data provided in 
IWalker et all (|2009l ) to remove potential interlopers with 
probability of membership less than 0.95 one retains a sam- 
ple of 2304 and 1351 stars for Fornax and Sculptor respec- 
tively. To ensure that there are at least two hundred stars 
in each bin which is a mini mum for a meaningful mea sure- 
ment of the fourth moment (|Amorisco fc Evan j[2012bl ') . the 
Fornax data is split into nine radial annuli of 230 stars plus 
an outer annulus of 234 stars and in the same manner the 
Sculptor data is split into five annuli of 225 plus one of 226. 
To construct the sampling distributions (see fig. 7 and sur- 



Table 1. Galaxy Parameters: Central coordinates (J2000), dis- 
tances and half-light radii 



Galaxy 




5c d [kpc] 


Ri/2 [kpc] 


Fornax 


02'"39™59'' 


-34°27.0' 138 


0.67 


Sculptor 




-33°42.5' 79 


0.26 



rounding text in [Richardson fc FairbairnI (|2012l ')') for Fornax 
and Sculptor (denoted with superscripts f and s respectively) 
we generated bootstrap sample sizes from a Pearson distri- 
bution of iV-'^ = 230 and iV= = 225 stars assuming scale dis- 
persions of cr/ = 1 0kms~^ and = 9kms~'^ and adding ex- 
perimental errors (|Amorisco fc Evansll2012"bl ) normally dis- 
tributed with widths 5^ = 0.22o-/ and 5° = 0.325o-f. Pro- 
jected stellar radii a r e calc ulated assuming galactic centres 
presented in iMateol (|l998l ) and converted to parsecs with 
distance measurements from that paper. The half-light ra- 
dius used in the Plummer model fit to the stellar density 
is taken from llrwin fc Hatzidimitrioul l|l995l ) and for clarity 
are shown in Table[T]with the aforementioned distances and 
central coordinates. 

4.2 MCMC Results 

4-2.1 Fornax 

Considering first the dispersion-only analysisin the upper 
panel of Fig. |4] we see that as expected the anisotropy is 
almost completely obscured by the traditional degeneracy 
of the flat dispersion measurement with the 95% confidence 
intervals at best enclosing j3{r) between moderately radial 
and strongly tangential orbits and placing even less con- 
straint at the galactic centre and larger distances. The mass 
contraints for the dispersion-only analysis are similar to the 
displayed error bars for the dispersion-only analysis of the 
same data in i Walker et al ] (120091 whic h uses more spatial 
bins, a more flexible Zhao (|ZhaollT996l ) density profile and 
assumes a Gaussia n distribution for Sg in the likel ihood. As 
shown in fig. 7 of iRichardson fc Fairbairiil ^2Q1± the lat- 
ter assumption is good provided the kurtosis is not more 
than mildly leptokurtic as is the case in Fornax. The depro- 
jected radius rls at which the logarithmic density slope of 
the stellar density profil e 7» = — 3 and w here the mass is 
robust to the anisotropy ijWolf et al.ll2010l l. is related to the 
projected half-light radius via rls « 1.23-Ri/2 for a Plum- 
mer profile and we see the characteristic pinch in the mass 
at r* 3 = 822pc. The mass slope confidence intervals in the 
upper left panel show that the dispersion-only analysis is 
unable to distinguish between the cusped NFW profiles in 
black and cored profiles that preserve the constant density 
with r = 3 out to around 750pc, at high significance though 
the NFW profile is mildly favoured. 

Adding the kurtosis to the analysis in the bottom panels 
of Fig. |4] we note first the dramatic improvement in preci- 
sion of the anisotropy parameter /3 at the half-light radius 
where there is maximal information from the data in the 
kurtosis measurem ent which is sensitive to the anisotropy 
jLokas et al.ll2005l '). Uncertainty in the anisotropy at the cen- 
tre of the galaxy is large which hopefully reflects that with 
little reference to the data at very small radii the freedom in 
the anisotropy proflle removes any spurious effects from the 
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Figure 2. MCMC reconstruction of mock Fornax data: The posterior distribution of dispersion-only analyses is marked in dashed blue 
while the dispersion-kurtosis analysis is solid-green. For brevity we collect the anisotropy parameters by plotting distributions of the 
anisotropy at the projected half light radius denoted /S^ij. Vertical red lines and red dots show the input parameters pj^. 




Figure 3. Confidence intervals for mock Fornax data: Columns 1-3 show the median (solid blue line) and 67/95% confidence intervals 
(respectively green and yellow shading) for the logarithmic mass slope, mass profile and anisotropy parameter as calculated from the 
derived posterior distributions output by the MCMC. Dashed red lines show true input parameters and in column 1 NFW profiles (dotted 
black) are plotted with three different scale radii, Ts = (200, 600, 2000)pc, for reference. Column 4 is the fourth order anisotropy, relevant 
only to the dispersion-kurtosis analysis. 
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0.5 1.0 1.5 2.0 0.0 0.5 1.0 1.5 0.0 0.5 1.0 1.5 2.0 0.0 0.5 1.0 1.5 2.0 



r [kpc] r [kpc] r [kpc] r [kpc] 

Figure 4. Confidence intervals for Fornax: Key as per figure |3] but witiiout red dashed lines to indicate tlie input parameter s. In column 
1 the data point and error bar corresponds to the mass slope and error derived for Fornax in lWalker &: Penarrubial 1 I2OIII ) and for its 
radius we use the average of the two sub component ha lf light radii. In column 2 the data points show the confidence intervals derived 
from the dispersion-only Jeans analysis in IWalker et al] lj2009. ') of the same data sets used here for M(300pc), M(Ri/2) and M{Ri) where 
Rl is the projected radius of the last star. 



assumption of a Plummer profile that limits the central be- 
haviour of the dispersion and fourth moment. In the fourth 
column the anisotropy between fourth order moments /3', 
constrained only by the kurtosis measurement, shows the 
new fourth order degeneracy as less affecting than the tra- 
ditional degeneracy in (3 in the plot to the left. For Fornax 
however it is not possible to distinguish between radially or 
tangentially biased kurtosis measurements. This could ex- 
plain why only moderate improvement is made on the mass 
and mass slope constraints from the top to bottom panels in 
the two columns on the left in comparison to the simulated 
data. Indeed the primary effect on the Fornax MCMC of 
including the kurtosis is to squeeze the confidence intervals 
with little change to the best estimates. 

In summary, for Fornax cores and cusps are still not 
distinguished at high significance but whilst the NFW pro- 
file is still in excellent ag reement, cored profiles as ex tended 
as the one obtained in (|Walker fc Penarrubiall201 j ) are m 
slight tension with the data which favours a density pro- 
file with an inner slope somewhere between the two. This 
resul t is also contrary to predictions from the virial theo- 
rem jAmorisco et al. 



l|jardel fc Gebhardtl 



201^ ) and orbit superposition methods 



20121 ') though the latter only provides 



likelihoods for NFW and a small set of cored profiles. 

One interesting improvement in the joint analysis is in 
the precision of the scale radius measurement which shows 
as a pinch in the mass slope confidence intervals. For cored 
profiles in particular having two references to the data in 
the joint analysis one envisages that this helps to distinguish 
this natural turning point in the moment profiles from the 
anisotropy. As the kurtosis is effectively independent of the 
scale density the degeneracy between the scale density and 
radius inherent to the overall scaling of the dispersion is 
partially broken. 



4.2.2 Sculptor 

The dispersion-only analysis of Sculptor shownin the top 
panels of [S] also shows a large degeneracy in the anisotropy 
parameter but here with a clear bias towards tangential ve- 
locity anisotropies. Both NFW profiles (dotted black) and 
profiles with cores out to the half fight radius are both with- 
ing the 95% confidence interval. The 67% confidence interval 
slightly favours the stee per NFW profile, in tension w ith the 
data point obtained in IWalker fc Penarrubial (|201ll ) using 
their multiple stellar population method. 

Caution is urged in placing too much emphasis on this 
second order result due to our choosing only 6 radial bins in 
order to later directly compare our second and fourth order 
analyses (as mentioned earlier, in order to reliably fit the 
kurtosis we require bins containing a relatively large num- 
ber of stars). Firstly, this relatively small number of bins 
compared to studies in the literature that adopt as many 
as y/N bins for the variance is not an optimal balance of 
spatial and statistical precision. With little reference to the 
outermost stars this choice places a strong emphasis on the 
outermost data point which suggests a rising velocity disper- 
sion whilst including more bins, as in IWalker et all l|2009l ). 
indicates that the dispersion then falls again at larger radii 
making it more compatible with cored fits. As such we take 
this result with a grain of salt and place more emphasis 
on the 95% confidence intervals that show the degeneracy 
that has been found in the literature. This is the reason for 
any small difference between the results of our second order 
analysis and other in the literature. 

The inclusion of the kurtosis for sculptor is rather more 
dramatic than for Fornax. Here we see that though the 
anisotropy at the galactic centre is obscured the joint anal- 
ysis makes a clear prediction for radial anisotropy. This is 
partly mirrored by the anisotropy between fourth order mo- 
ments j3' though at large radii the significance is smaller 
owing to the fourth order degeneracy. More interesting still. 
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Figure 5. Confidence intervals for Sculptor: Key as per Fig. [4] 



the cored solution, in excellent agreement with the pre- 
diction from multiple stellar populations, is favoured over 
the NFW one to high significance up to and beyond the 
half-light radius breaking the degeneracy. This is also mani- 
fest in the mass estimate where systematically lower masses 
are predicted at small radii which would reduce the as- 
trophysical J-factor over small integration angles. On face 
value then we validate the prediction of a cored density pro- 
file in studies of mu ltiple populations (|Walker fc Penarrubial 
I2OIII: Amorisco fc Evan s 2012a ) and from the virial theorem 
(|Agnello fc Evansll201 j r. 



with a similar DM scale radius logj,, r_ 



2.81 



+0.17(1.28) 
-0.07(0.13) 



4.3 Why is the Fourth Order Method Better for 
Sculptor and not for Fornax? 

The results for Sculptor and the simulated data imply that 
the kurtosis measurements are key to distinguishing between 
the scale radius r_2 and shape parameter a, a degeneracy 
that normally exists with a fiat velocity dispersion profile, 
however there is no corresponding significant improvement 
for Fornax. Inspection of the kurtosis data in Fig. [S] shows 
that while Fornax has an approximately flat series of kur- 
tosis measurements that dips below the Gaussian value of 
AC = 3, the sculptor data exhibits a group of leptokurtic 
points (k > 3) in the dense stellar region before dipping 
again to Gaussianity at t he e nd confirming the res ults in 
lAmorisco fc Evand l|2012bl ) and lBreddels et al.l (|2012l ). With 
the assumption of mild anisotropy we see that the kurtosis 
profiles converge at a point below the half-light radius. The 
location of this convergence point on the kurtosis axes dis- 
plays a dependence on the shape parameter a, an increase 
in which shifts it upwards. 

Comparing solutions with identical DM density profiles 
in the upper and lower panels of Fig.[S]we see that by extend- 
ing the half-light radius in the Jeans equation the scale of 
the shift is reduced and the convergence points are less sep- 
arated on the kurtosis axis making it harder to distinguish 
between them. For Fornax where the median DM scale ra- 
dius returned by the MCMC logjo r-2 = ^.Sl^gJgjp'jQj is 
similar to the projected half-light radius r_2/i?i/2 « 1 we 
see that the shift is reduced compared to Sculptor which 



has a more embedded stellar population r_2//?i/2 ~ 2.5. 
The NFW-like density profile with a = 0.2 is least sensitive 
to this effect and a very large scale radius(and thus diffuse 
halo to fit the dispersion data) is required to shift the con- 
vergence point towards the Gaussian value k = 3. For this 
reason it is very difficult to generate leptokurtic points such 
as those in the Sculptor data set with k > 3.5 that lie near 
the convergence radius. This is exacerbated by the fact that 
fat tailed platykurtic distributions are less prone to large sta- 
tistical fiuctuations which reduces the likelihood of points in 
the outer tails. 

In conclusion there appears to be a radius where the 
kurtosis is invariant to the anisotropy which for NFW-like 
profiles falls well below k = 3. While radial anisotropy shifts 
the curve upwards towards the Sculptor data points at larger 
radii, the kurtosis between 100 and 150pc is stuck and can- 
not describe the data points with k > 3.5. For this reason 
not only is radial anisotropy preferred in Sculptor but the 
MCMC is able to fit cored profiles with higher likelihood 
than steeper ones regardless of the anisotropy. In Fornax 
the data points have kurtosis k < 3 and are easier to de- 
scribe with NFW-like profiles. The extended half-light ra- 
dius also reduces the ability to distinguish between DM den- 
sity profiles at the point where the kurtosis is invariant to 
the anisotropy. If these empirical results stand up to a wider 
range choices for the anisotropy parameters or a variation 
in the adopted form of the stellar density profile then Sculp- 
tor's leptokurtic data points could be a smoking gun for a 
cored density profile and we are motivated (Richardson fc 
Fairbairn in preparation) to confirm these findings with an- 
alytic reasoning. 



4.4 Key Assumptions and Limitations 

In this subsection we address the potential impact of the 
simplifying assumptions inherent to the Jeans analysis used 
above as well as astrophysical processes and phenomena 
that, for simplicity, were omitted. 

The assumption of spherical symmetry, necessary to 
close the set of moment equations, has been shown 
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Figure 6. Generic kurtosis profiles and data: Kurtosis l|26|l is 
plotted for the ten radial bins of Fornax (lower panel) and 6 of 
Sculptor (upper panel). Against these we plot line of sight kur- 
tosis profiles from the Jeans equations assuming Plummer pro- 
files for each with half-light radii from Table 1 and Einasto scale 
radii of r_2 =700pc similar to the median scale radii r_2 =698 
and 641pc returned by the fourth order MCMC for Fornax and 
Sculptor respectively. With no attempt made at fitting the data 
we show three constant anisotropy models with /3 = (—0.6, 0, 0.4) 
and isotropy at fourth order 0' = for each adopted shape pa- 
rameter a that is distinguished by colour and line style. A simple 
means to distinguish between anisotropy is to note that tangential 
curves diverge at the centre, isotropic curves tend to a finite con- 
stant value and radially anisotropic curves vanish at the origin. 
Curves in both plots have the same DM density and anisotropy 
parameters with differences arising because of the change in the 
adopted half-light radius. 



l|Evans et al.ll2009l ) to have a significant impact on the dy- 
namics at the Galactic centre when coupled with restrictions 
on the stellar and dark matter logarithmic density slopes. 
Though concern is raised that this artifact of the assump- 
tions could lead to spurious results as discussed in Section 
3.1 we take comfort from the fact that the uncertainty in the 
anisotropy parameter derived from the MCMC is large at the 
centre, thus minimising its influence with the anisotropy in 
the densest stellar regions free to vary independently from 
these constraints. For the joint analysis, we note that the 
confidence intervals narrow considerably at a finite radius 
near that of half light where there are most stars and thus 
information from the data and that for Sculptor, with the 
strongest constraints this effect is particularly pronounced. 
This suggests that adopting a generalised radial profile for 
the anisotropy is important to decouple this useful behaviour 
from the assumptions though to confirm these findings one 



should perform an analysis where both inner slopes have a 

greater degree of freedom. 

Photometric measurements (|lrwin fc Hatzidimitrioul 

I1995I ) of Fornax and Sculptor indicate a departure from 
spherical symmetry with the ratio of major and minor axes 
~ 0.7. By consid ering elliptical rather than ci rcular annuli it 
has been shown l|Walker fc Penarrubia|[2"oill ) that the mass 
slopes derived by studying multiple populations, where un- 
certainties in the half light radii are important, are sensitive 
with Sculptor dropping f rom F = 2.95 to F = 2.4. By con- 
trast Amorisco fc EvansI l|2012bl ) find that binned measure- 
ments of the sample moments are broadly similar in the mi- 
nor and major axis directions and with more spatial points 
of reference the Jeans analysis should be more robust. 

The other major assumption inherent to the Jeans anal- 
ysis is dynamic equilibrium though significant tidal dis- 
ruption in the classical dwarfs h as only been conclusively 
identified jMaiewski et al. 20031) in Sagittarius. Though 



simulations (iKlimento wski et al.l |2007| ) suggested that the 
Jeans analysis is robust to tidal forces with suitable inter- 
loper removal procedures, it has recently been suggested 
that measurements obtaine d with the mass slope method 
(| Walker fc Penarrubiall201ll ') are subject to systematic un- 
certainties dependent on the unknown line of sight. One of 
the great powers of this method, however, is that it can be 
applied to small samples and an analysis of a larger group 
of data sets would reduce any systematic bias. 

Bin ary stars cou ld also inflate the dispersions 
( Hargreaves et aLlll996l) though th e case has been put for- 
ward ( McConnachie fc Cotel |201G| ) that significant effects 
only occur for dispersions on the scale of the ultra faint 
dwarfs with ap < 3kms~'^ significantly smaller than those of 
the classical dSphs considered herein with ap ~ lOkms"^. 



5 DISCUSSION 



In iRichardson fc FairbairnI (|2012| ) a Jeans analysis of the 
variance and kurtosis was devised that for the first time 
fully extended the framework from second order to fourth 
order without imposing restrictions on the anisotropy be- 
tween fourth order moments which we introduced as /3'. 

With generic and uncorrelated parametric forms for the 
radial profiles of /3 and P' we applied the analysis to kine- 
matic data sets for the Milky Way dwarf spheroidals Fornax 
and Sculptor in the hope of breaking the mass-anisotropy de- 
generacy that plagues dispersion-only analyses of real data 
sets with flat dispersion curves. The assumption of spheri- 
cal symmetry and a fixed parametrised form for the stellar 
density profile iy{r) has a significant impact upon the Jeans 
equations at the galactic centre and the freedom to divorce 
the central anisotropy from that at regions of high stellar 
densi ty provides more confidence than previous joint anal- 
yses (|Lokas et al.ll2005h that have been limited to constant 
anisotropy. 

Though the joint analysis of Fornax was able to sub- 
stantially improve the precision in the anisotropy measure- 
ment, little improvement was made for the mass slope which 
could not distinguish between an NFW-like Einasto profile 
and one with an extended core at high significance. Slight 
tension was however found with large cores predicted in the 
literature from studies of multiple populations with the best 
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fits being intermediates between these and the NFW-hke 
models. By contrast the analysis of Sculptor favoured radial 
anisotropy and showed a clear preference for cored solutions 
in excellent agreement with the literature. At first glance it 
appears that the rather distinct kurtosis measurements for 
the two galaxies are key to interpreting these results and 
we demonstrate that the leptokurtic data points in Sculptor 
seem difficult to reproduce with an NFW-like profile when 
one assumes moderate anisotropy.Purther investigation into 
this effect with a wider range of anisotropy and stellar den- 
sity models is required to confirm this finding which we hope 
to present in the near future. As both galaxies display two 
very different anisotropy and density profiles it could be ar- 
gued that astrophysical effects, more chaotic in nature, play 
a significant role in shaping the density profiles of dwarf 
spheroidal galaxies that could obscure systematic effects in- 
duced by warm or self interacting dark matter cosmologies 
for which a larger sample of galaxies is needed. 

In summary this method provides an analytical insight 
into the dynamics of spherical systems with a generality bet- 
tered only by numerical orbital superposition methods. It 
also complements the bold predictions from stellar sub com- 
ponent methods by testing the assumption of Gaussianity. 
With further investigation into the nature of the kurtosis 
profile we hope to provide simple analytic insights that, like 
the robust mass estimate at the half-light radii, could give 
simple justification to the interesting results of the MCMC 
analysis presented here and help to settle the cusp vs core 
debate. 
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